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ABSTRACT 

Formulas  for  computing  the  gravitational  effect  of  some  simple  two- 
and  three-dimensional  geometric  figures  are  presented  in  forms  suitable 
for  use  with  digital  computers  and  in  many  cases,  programmable  desk 
calculators.  Basic  con5)utation  schemes  are  presented  for  con^jlex 
two-  and  three-dimensional  bodies  of  arbitrary  shape.  Some  simple 
inversion  rules  or  techniques  are  presented  which  yield  approximations 
of  depth  based  on  simple  geometric  figiires.  Such  inversion  techniques 
are  particularly  applicable  where  gravity  measurements  and/or  other 
geophysical  data  are  sparse.  These  techniques  yield  first  approximations 
of  the  depth,  size,  and  shape  of  the  mass  or  masses  causing  a given 
residual  gravity  anomaly.  Density  or  density  contrast  is  independent 
of  depth,  size  and  shape  and  is  discussed  separately.  Density  and 
porosity  are  defined  by  appropriate  equations.  Generalized  density 
relationships,  based  on  rock  types  are  discussed.  Equations  for 
determination  of  effective  density  from  gravity  measurements  are 
presented.  The  concluding  remarks  cover  some  general  aspects  of 
gravitational  modeling. 
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3-1  Two-Dimensionel  Criterion 


LIST  OF  NOTATIONS 


The  notation  is  uniform  insofar  as  possible.  Any  change  in  the  notation 
listed  below  is  specifically  noted.  The  notation  is  generally  based 
on  the  following  definitions: 


X,  Y,  Z 


Ax,  Az 


= Define  the  three  coordinate  axes  of  a left-handed 
Cartesian  coordinate  system  in  three-dimensional 
models . 


= Define  the  two  coordinate  axes  in  two-dimensional 
models. 


= The  horizontal  ground  distance  along  the  X-axis, 
usually  from  a point  above  the  center  of  the  body 
to  the  computation  point. 


= The  horizontal  ground  distance  along  the  Y-axis, 
usually  from  a point  above  the  center  of  the 
body  to  the  con^JUtation  point. 


The  depth  from  the  surface , usually  the  X , Y plane , 
to  the  center  of  the  body. 


= The  distance  from  the  center  of  the  body  to  the 
computation  point. 


= The  radius  of  the  body. 


= Thickness  of  laminar  bodies. 


= Universal  constant  of  gravitation  = 6.673  x 10“® 
cm®/(grams) (sec)^ . 


= Calculated  gravitational  effect,  usually  the  vertical 
component  of  the  acceleration  of  gravity  (in  milligals) 


Gravitational  potential 


= Volime  density  or  density  contrast  as  a function 
of  X,  y,  z. 


= Siirface  density  or  density  contrast  as  a function 
of  X,  z. 


= Cross-sectional  area  of  two-dimensional  bodies  as 
a fvinction  of  x,  z. 


Volume . 


» Mass. 
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1.  INTRODUCTION 


The  ptirpose  of  this  publication  is  to  present  a unified  discussion 


of  several  widely  accepted  analytical  techniques  of  gravitational 


modeling  along  with  some  elementary  techniques  of  interpretation  in 


the  form  of  a reference  manual.  The  formulation  of  these  techniques 


is  oriented  towards  applications  compatible  with  programmable  desk 


calculators  and  digital  computers. 


Gravitational  modeling  refers  to  the  direct  problem  o*  determining 


the  mass  attraction  of  a geologic  body  or  structure  of  known  t,hape. 


depth,  size,  and  density  or  density  contrast.  The  reliability  of  the 


solution  depends  on  how  well  the  shape  of  th'!  body  or  structure  is 


known  and  how  closely  it  can  be  approximated  by  analytical  expressions. 


Thus,  the  direct  problem  can  be  solved  only  when  specific  relatively 


homogeneous  geologic  bodies  or  structures  can  be  identified  as  the 


gravitating  masses. 


Such  detailed  knowledge  of  geologic  structure  is  limited  to  features 


in  the  upper  portion  of  the  earth's  crust.  The  attraction  of  such 


bodies  is  assumed  to  be  an  approximation  of  the  local  or  residual  com- 


ponent of  the  observed  gravity  anomaly  field. 


Some  common  structiires  which  significantly  influence  the  residual 


gravity  anomaly  field  are:  well  defined  topographic  features,  sedi- 


mentary basins,  faults,  anticlines,  synclines  and  intrusive  bodies  of 
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discernible  shape. 


Gravimetric  interpretation  refers  to  the  inverse  gravimetric  problem. 
The  solution  to  the  inverse  problem  has  as  its  objective  the  deter- 
mination of  the  unknown  parameters  (i.e.  , shape,  depth,  and  density 
or  density  contrast)  from  a known  residual  gravity  anomaly  field.  A 
unique  solution  to  the  inverse  problem  is  theoretically  impossible 
because  a given  residual  gravity  anomaly  field  can  be  produced  by  an 
infinite  nxamber  of  mass  conf igxarations . 

However,  reliable  first  approximations  of  the  unknown  parameters 
can  often  be  obtained  from  the  shape  and  magnitude  of  residual  gravity 
anomaly  profiles  as  well  as  from  various  other  sources  of  geologic  dnd 
geophysical  data.  Such  sources  include:  field  maps,  electric  logs, 
seismic  profiles,  and  borehole  information.  The  first  approximations 
of  the  unknown  parameters  can  then  be  used  in  an  iterative  inversion 
procedxire.  In  such  a procedure,  the  unknown  parameters  are  then  suc- 
cessively adjusted  until  an  observed  residual  gravity  anomaly  field 
is  approximated  within  some  predetermined  tolerances.  The  adjusted 
unknown  parameters  must  then  be  examined  within  the  context  of  geologic 
and  geophysical  realism.  It  is  often  necessary  to  fiarther  constrain 
the  allowable  range  of  any  or  all  of  the  unknown  parameters  in  suc- 
cessive models. 

The  inversion  of  residual  gravity  anomalies  can  yield  no  better 
results  than  the  residual  gravity  anomalies  themselves.  There  are 
certain  inherent  ambiguities  in  the  separation  techniques  for  computing 
residual  gravity  anomalies.  The  ambiguities  are  due  to  the  fact  that 


f 
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measioring  techniques  depict  the  influences  of  all  the  masses  within 
the  measuring  range  of  the  instr'.unents.  Thus,  it  is  impossible  to 
completely  separate  the  regional  and  residual  anomaly  fields.  Therefore, 
care  must  be  taken  to  minimize  observational  errors  and  provide  the 
best  possible  observed  data. 

This  report  contains  a relatively  complete  list  of  the  formulas 
and  computing  schemes  which  form  the  basis  for  the  solution  of  the 
direct  problem.  However,  there  is  only  a partial  list  of  formulas 
for  depth  and  density  determination.  Such  fcrmv'las  yield  realistic; 
ranges  of  depth  ana  densi  ty  and  do  P'-i.  conFti  e a solution  of  the 
inverse  problem.  Hov-n'ci*,  they  can  ..'ic'ld  valid  first  approximaticjns 
of  depth  and  density  for  use  in  gravii  ational  modeling.  Detailed 
approaches  to  the  solution  of  the  inverse  problem  are  be/ond  the  Tjuri^osc 
and  scope  of  this  pubLi  cation  and  can  be  found  in  many  <>!'  tur  ci  ted 
references . 

The  two-dimensional,  and.  threo-dii!.f.:,Tsi':nal  fomiuias  are  pies'-ntcci 
in  the  general  formats  used  by  Heiland  (1968)  and  Taiwan!  (1971).  Tin. 
functional  notation  usc:d  icith  the  three-dimensional  formulas  is  identi- 
cal with  that  presented  by  talwar.i  (19T^). 

1 . 1 Mathematical  Development 

1.1.1  The  ^iree-Dimonr.ic.nal  Cas'' 

Based  on  Newton's  inverse  souorc  law,  the  gravitational 
potential  U of  a gravitating  bodj"  of  Volume  V (Figure  1-1 ) is  given 
by: 
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U = K p 


(1-1) 


where : 


r - {x2  + y2  + z^Vh 

dV  = dx  dy  dz  (volxune  element.  Figure  1-2) 


(1-2) 


Expressed  in  rectangular  coordinates,  the  gravitationea 


potential  of  V becomes: 


U = K p 


dx  dy  dz 
(x^  + y^  + z^J 


(1-3) 


In  spherical  coordinates,  r,  <p,  a (Figvire  1-3),  the  gravi- 


tational potential  is  then  expressed  by: 


U = KpjJJrdr  cos(J>  dcj)  da 


(1-U) 


Equation  1-4  is  derived  by  malting  the  following  substitutions 

in  equation  1-3: 

r = (x^  + y2  + z=^)i/2 

r^  cos<J)  dr  d(|)  da  = dx  dy  dz  (volume  element  (Figure  1-2)) 
r sin({)  = z 

The  gravitational  attraction,  g^,  of  a gravitating  body 
of  volume  V is  the  partial  derivative  of  the  potential  U with  respect 
to  z,  and  is  defined  as  follows: 


Sv  = T7  = ^ P 


(1-6) 


t 


Figures  1-2 

Volume  Elements 

Cartesian  Coordinates 


dR' 


Cylindrical  Coordinates : 


Figure  1-3 
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Coordinate  Systems 


Expressed  in  rectang\ilar  coordinates.  Figure  1-1,  equation 


1-6  becomes; 


a , p j J j V dx  dy  dz  (1-7) 


In  a spherical  coordinate  system,  r,  (j),  a (Figure  1-3), 
the  gravitational  attraction  is  given  by: 


g^  = K p I I sin(J)  cos({>  dr  d()>  da 


(1-8) 


Equation  1-8  is  derived  from  eq\iation  1-T  by  making  the 
following  substitutions: 

r>  = (x“  t y"  ♦ 

r*  cos<|)  dr  d<|>  da  = dx  dy  dz  (Figure  1-2)  (l-9) 

r sin<t>  = z 


Expressed  in  a cylindrical  coordinate  system,  R , z,  ot 
Figure  1-3,  the  gravitational  attraction  is  given  by  equation  1-10: 


„ - r f f R'  z dz  dR^ 

P JJ  J 


da 


(1-10) 


where,  equation  1-10  is  derived  from  equation  1-7  by  the  following 
substitutions . 

(R'*  + z^)V2  = + y^  + z^)V2 

R'‘z  dz  dR'’da=  z dxdydz  (Figiure  1-2) 

In  a conical  coordinate  system  R'’,  <|)  ,a  (Figure  1-3),  the 


(1-11) 
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gravitational  attraction  is  then  given  by  equation  1-12. 


g = K p 

V 


sin(J)  dR  " d({i  da 


(1-12) 


Equation  1-12  is  derived  from  equation  1-7  by  madting  the 
following  substitutions; 

= (x^  + y2  + z2)  (1-13) 

dR'  d(|>  da  = dx  dy  dz 
R''  sin4»  = z 


1.1.2  The  Two-Dimensional  Case 

Many  topograph:' c md  geologic  rlructures  are  elongated  along 
their  strike.  Such  structures  can  often  tr-  r^od.eled  by  two-dimensional 
approximations.  The  following  paragraphs  present  the  basic  mathematical 
development  of  the  two-dimensional  approach. 

Two-dimensional  structures  are  modeled  by  positioning  the 
coordinate  axes  such  that  the  Y-axis  is  infinite  in  extent  and  is 
parallel  to  the  strike  length,  or  elongated  dimension,  of  the  structure 
to  be  modeled.  The  shape  of  the  structure  is  then  defined  by  a vertical 
cross-section  of  area,  A^,  Figure  1-U. 

The  formula  for  computing  the  gravitational  attraction  of 
two-dimensional  bodies  is  derived  by  integrating  equation  1-T  with 
limits  from  -«>  to  +<*>  in  Y. 


g^  = Kp 


z dx  dz 


+03  dy 

-oo  (x^  + y^+ 


(l-l4) 


The  integral  to  be  evaluated  is  then  rearranged  and  becomes; 


9 


i 


i 


BEST  AVAILABLE  COPY 


i I 


_ f"*’'*’ 


J2_ 
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(x2  + (1  ^ ■ y-  ■;) 


The  following  substitutions  are  made; 


(1-15) 


tan  u = / -,'5 177  and  = (x^  +12)1'/!^ 


(1-16) 


Then,  equation  1-15  becomes: 

+Tr/2 


Tr/2 


cosu  du 


(1-17) 


■ (x^  + zM 

Equation  1-17  is  substituted  into  equation  1-lU  and  the 


resultant  formula  is: 


K p f 

J A 


TIFT-PT  ^ 


(1-18) 


Equation  I-I8  expressed  in  polar  coordinates  r,  ([)  becomes, 
(Figure  1-1+) : 


g _ = 2 K p I I sincj)  dr  dij) 


(1-19) 


Two-dimensional  formulas  can  be  made  to  closely  approximate 
bodies  of  finite  extent  by  making  "end  corrections"  for  the  finite 
strike  length,  y (Figure  1-1+ ).  The  "end  correction"  for  an  elongated 
body  of  finite  strike  length  is  defined  by  the  following  integral 
equation: 

ft  z dx  dz 

(1-20) 


•Sg,,  = 2 K p y 


(x*  + z^)  /x""  + z^  + y^ 
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which  is  given  in  polar  coordinates  by: 


, , K p y I f (1-21) 

^ ■'A-'  + 

Equations  1-20  and  1-21  cannot  be  integrated  in  closed  form 
and  mvist  be  evaluated  numerically.  The  vertical  component  of  the 
attraction  (i.e.,  gravitational  attraction)  of  the  bodies  disciissed 
in  following  sections  will  be  referred  to  as  "gravitational  effect." 


( 


2.  THREE-DIMENSIONAL  ATTRACTION  FORMULAS 
2. 1 Some  Simple  Three-Dimensional  Formulas 

All  the  formulas  presented  in  this  section  eire  derived  from  equa- 
tion 1-6  of  the  proceeding  section  on  mathematical  development. 

Talweini  (1973)  shows  that  majiy  of  the  formulas  can  be  expressed  in 


functional  form  as  follows: 


= K p f(x,y,z)  F(x/z) 


(2-1) 


where  f(x,y,z)  is  a variable  function  of  x,  y,  or  z and  F(x/z)  is  a 
dimensionless  scaling  function.  The  function  F(x/z)  is  dependent  on 
the  variation  of  x/z  from  zero  to  one.  I'-erefore , the  evaluation  of 
the  formulas  can  be  greatl;,'  ^iTiplified  u;.;  ■ riots  of  F(x/z)  against 
the  ratios  x/z  and  z/x. 

2.1.1  Rectangular  Volume  Element 

The  dimensions  of  the  volume  elcr.  .rnt  are  given  as  Ax,  Ay, 
and  Az  as  shown  in  Figure  2-1.  The  gravitational  effect  of  the  volume 
element  is  given  by  equation  2-2. 


g^  = K p Ax  Ay  Az  ^ 


where : 


= K p Ax  Ay  Az  z 


r = (x^  + y^)'/^ 


(r^  + z^)  = 


Tho  function  Fi(r/z)  is  then  defined  as: 

Fi(r/z)  = (1  + 

and  is  plotted  in  Figure  2-2. 


BtSI  AVAIUBLE  COPY 


(2-2) 


(2-3) 


The  volume  element  is  not  very  useful  in  modeling  real 
geologic  bodies  or  structures.  It  is  included  in  this  discussion  for 
the  sake  of  completeness. 

2.1.2  Sphere 

The  gravitational  effect  of  a solid  sphere  is  derived  by 
substituting  the  mass  of  a sphere  for  the  mass  of  a volume  element  in 
equation  2-1.  Equation  2-k  gives  the  gravitational  effect  of  a solid 
sphere  as  shown  in  Figure  2-3. 


= “I  TT  K p z (R/r)^ 


(21U) 


= ^ IT  K p (rVz^)  Fi(ri/z) 


where ; 


•i  = (x^  + y^) 


(2-5) 


The  fimction  F^(r^/z)  is  given  in  Figure  2-2. 

The  mass  of  a solid  homogeneous  sphere  is  distributed 
symmetrically  about  the  center.  Thus,  the  total  and  horizontal  components 
of  the  gravitational  attraction  are  easily  computed  as  shown  in  equations 
2-6  and  2-7. 


g^  = Total 
&J,  = gy  (r/z) 


(2-6) 


1 


g = Horizontal  Component 

XI 


gR  = 6r  (ri/r)  (2-7) 

= 6p  sine 

The  equation  for  the  gravitational  effect  of  a sphere  is 
identical  with  that  of  an  equivalent  point  mass  concentrated  at  the 
center  of  the  sphere.  The  spherical  model  can  be  used  to  juDdel  masses 
of  simple  geometry  or  masses  in  which  the  geometry  is  poorly  defined. 
Equation  2-k  is  the  basis  of  depth  range  determination  as  will  be  dis- 
cussed later.  Equation  2-h  also  forms  the  basis  of  the  so-called 
"point  mass"  approach  to  gravity  field  modeling. 

2.1.3  Horizontal  Line  Element 

The  horizontal  line  element  with  a cross-sectional  area  ^ 
is  parallel  to  the  Y-axis  and  extends  from  yi  to  y2  as  shown  in  Figure 
2-h.  Equation  2-8  then  gives  the  gravitational  effect. 


M 

g = K p — 
V z 

(1  + xVz2)“l 

172  Yi' 

(2-8)  r 

= K 

fyz 

Fz  (x/z)  — 

i 

z 

\r2 

I 

j 

where ; 


ri  = (x^  + yi^  + z^ ) 

r2  = (x^  + y2^  + z^)*/^  (2-9) 

F2(x/z)  = (l  + x^/z^)“* 

The  function  F2(x/z)  is  plotted  in  Figure  2-5. 
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The  problem  becomes  a two-dimensional  problem  when  the 
line  element  is  infinite  in  length  (i.e.,  extending  from  y = -<»  to 


y = +“).  Equation  2-8  is  then  reduced  to: 


g^  = 2 K p ^ F2(x/z) 

V Z 


(2-10) 


Equation  2-10  also  applies  to  an  infinite  horizontal  cylinder  of  cross- 
sectional  area  AA. 

The  horizontal  line  element  is  used  mainly  to  determine  the 
two-dimensional  criterion  as  will  be  shown  in' later  sections. 

2.1.1*  Vertical  Line  Element 

The  vertical  line  element  of  cross-sectional  area  M is 
oriented  parallel  to  the  z-axis  and  extends  from  zi  to  zj  as  shown  in 
Figure  2-6.  The  gravitational  effect  is  given  by  equation  2-11. 


g^  = K p M 


= K p AA 


1 1 

(x^  + y^  + z 1^ ) ~ (x^  + y^  + Z2^ ) ^ ^ 


(r^  + zi^)  (r^  + Z2^ ) 


(2-11) 


_ K (j.Qg0j  _ COS02) 


where : 


r = (x^  + y^)‘/^ 
zi 
r 
Z2 
r 


01  = teui"  — 

02  = tan"^ 

If  Zi  = 0,  then  equation  2-11  simplifies  to: 


(2-12) 
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= K p M 


'l  

r ~ (r^  + z: 


(2-13) 


K p M ,,  „ , 

= (1  - cos  62) 


If  Z2  = ",  then  equation  2-11  becomes: 
= K p M [(1.2  + 


= K p M 


If  zi  = 0 and  Z2  = ",  equation  2-11  is  further  simplified  to: 


(2-IU) 


K p M 


(2-15) 


The  equations  for  the  gravitational  effect  of  a vertical  line  element 
can  be  used  to  approximate  the  attraction  of  right  vertical  prisms 
and  cylinders  (i.e.,  volcanic  plugs  and  salt  domes),  as  long  as  the 
cross-sectional  areas,  M,  are  small  compared  to  the  other  dimensions 
of  the  featvire. 

2.1.5  Vertical  Rectangular  Lamina 

In  Figiire  2-7,  the  vertical  lamina  ab  is  defined  by  the 
opposite  corners  a(x,0,0)  and  b(x,y,zi).  Lamina  ab  is  oriented  peurallel 
to  the  YZ  plane  with  dimensions  y,  zi  and  Ax.  The  gravitational  effect 
of  lamina  ab  is  given  by  equation  2-l6. 


Sv  ' ^ P T ^ 


[{(x^  + y^  + zi^)‘/^  - y}  {(x^  + y^)^/^  + y}" 


= K p Ax  In  ' ■ 


(x^  + y"  + zi^)*/"  + y}  {(x^  + y^ 


r(x^  + zi)*/^  (y  + (x^  + y^)^'^^}' 


(2-16) 


X ly  + (x  + y + zi‘ 


The  natural  logarithm  function  can  be  expressed  as  follows: 
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Figlire  2-8 


r(l  + ZiVx2)l/2  {1  + (1  + x2/y2)l/2}] 

G(y/x,Zi/x)  = In  1 + (l  V x'^'/y'S' T z ^ Vy^')'^'/^ 


then: 


= K p Ax  G(y/x,zi/x) 


(2-18) 


The  function  G(y/x,Zi/x)  is  plotted  in  Figure  2-8. 

The  gravitational  effect  of  an  arbitrary  rectangular  lamina 


gd  with  a top  edge  at  z^  and  bottom  edge  at  Z2  is  given  as: 
g^  = K p Ax  [G(y/x,Z2/x)  - G(y/x,Zi/x)] 


(2-19) 


For  a vertical  lamina  of  infinite  extent  in  y and  ranging 
in  vertical  extent  from  z = 0 to  z = Zj,  equation  2-16  is  simplified 


|x^  + 


g = K p Ax  In 

V 


=KpAxln(l+Zi^ /x^ ) 


(2-20) 


= K p Ax  F3(zi/x) 


where : 


F3(zi/x)  = In  (l  + zi^/x^) 


(2-21) 


The  function  F3(z/x)  is  plotted  in  Figure  2-9. 

If  the  infinite  vertical  lamina  ranges  in  vertical  extent 
from  z = zi  to  z = Z2,  the  gravitational  effect  is  given  as: 


g^  = K p Ax  [F3(z2/x)  - F3(zi/x)] 


= K p Ax  In 


X^  + Z2^ 

n2  .■  Z 
X + z 1 


(2-22) 


The  vertical  rectangular  lamina  is  best  used  to  model  vertical  or  near 
vertical  dikes.  Functions  G and  F3  are  unsuitable  for  small  values  of  x. 
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Figure  2-9 
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Figure  2-10 


Horizontal  Rectangular  Lamina 


2.1.6  Horizontal  Rectangular  Lamina 

The  horizontal  lamina  ah  in  Figure  2-10  is  oriented  parallel 
to  the  XY  plane  with  dimensions  x,  y,  and  Az.  Equation  2-23  then  gives 
the  gravitational  effect  of  lamina  ah. 


= K p Az 


IT  . 

- - sin  1 


yi 


- sin 


- 1 


(x^  + z^ ) (x^  + Y ^ 

(yi^'+  z2)i/2  -(72  /z-zyr/T 


(2-23) 


= K p Az 


YiZ 


2 X (x^  + y 1^  yi(x^  + yi^  + 


xz 


= K p Az 


tan' 


_i 


xyi 


(x‘‘  + y 


Then  the  function  in  brackets  in  equation  2-23  can  be  expressed  as 
follows : 


xyi 


G^Cxyi/xyiz)  = tan"!  ^ ^ ^ \ 


(2-21+) 


Therefore  the  gravitational  effect  of  lamina  ah  is  given  as : 


g^  = K p Az  G2(xyi/xyiz) 


(2-25) 


The  gravitational  effect  of  an  arbitrary  horizontal  lamina  such  as  cf 
is  then  given  as: 


g^  = K p Az  [G2(xy2/xy2z)  - G2(xyi/xyiz)  + G2 (xiyj/xjyiz)  - G2(xiy2/xiy2z) ] 

...(2-26) 


where  G2'’  = (tt/2)  - G2  is  plotted  in  Figure  2-11.  The  horizontal  rec- 
tangular lamina  is  used  to  approximate  sills  or  layered  sedimentary 
features . 
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Figure  2-11 


Function  Q^'  (x^y^/x^y^z) 
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Figure  2-13 
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Thin  Circular  Horizontal  Disk, 
General  Expression 


2.1. T Thin  Circular  Horizontal  Disk,  Points  on  the  Axis  ThroTigh 


the  Center 


The  disk  is  defined  in  a cylindrical  coordinate  system  as  shown 

in  Figure  2-12.  The  gravitational  effect  of  the  disk  is  given  by  equation 

2-27.  j- 

g^  = 2TiKpA2  1- 

= 2 IT  K p Az  (l-sin0) 

2.1.8  Thin  Circular  Horizontal  Disk,  General  Expression 

Again  the  disk  is  defined  in  a cylindrical  coordinate  system, 

o 

see  Figure  2-13.  The  gravitational  effect  of  the  disk  is  given  by  equation 


2-28. 


= K p Az  0)  (2-28) 

The  parameter  oj  is  defined  as  the  solid  angle  subtended  by  the  median  plane 
of  the  disk  at  the  computation  point  P.  The  solid  angle  co  has  been  deter- 
mined empirically  and  published  in  template  form  by  a number  of  investiga- 
tors: Nettleton  (l9Tl)  and  Talwani  (1973).  The  solid  angle  template  in 
Figure  2-lit  is  reproduced  from  Talwani  (1973)  with  the  permission  of  the 
author  and  publisher.  The  template  gives  solid  angles  for  solid  circiilar 
disks  and  cylinders  as  a fxmction  of  the  ratios  x/R  and  R/z. 

Equation  2-27  is  an  approximation  based  on  the  assumption  that 
al 1 mass  of  the  disk  or  cylinder  is  concentrated  on  the  median  plane. 
Nettleton  (197 l)  shows  the  approximation  to  be  in  error  by  -2%  for 
Az  < (l/2)  z and  ^Iz  = 1. 

The  thin  disk  has  the  same  application  as  the  horizontal  rectangular 
lamina  with  the  added  advantage  of  computational  ease. 
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Figure  2-lU 


Solid  Angle  Chart  - Function  Gj^ 
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2.1.9  Vertical  Right  Circuleo'  Cylinder  of  Finite  Depth , Point  on 
the  Axis 

The  right  circtilar  cylinder  in  Figure  2-15  is  defined  in  a 
cylindrical  coordinate  system.  Equation  2-29  then  gives  the  gravitational 
effect  of  the  cylinder: 

g^  = 2 TT  K p [(Z2  - Zi)  - (Z2^  + + (zi^  + R^)!/^]  (2-29) 

= 2iTKp{h-a  + b) 

where : 

h = (Z2  - Zi) 

a = (z22  + R2)i/2  (2-30) 

b = (zi"  + r2)V2 

If  Zj  = 0,  then  equation  2-29  is  simplified  to: 

g^  = 2 IT  K p [R  + Z2  - (R^  + Z2^)V^]  (2-31) 

For  Z2  = equation  2-29  then  becomes: 

g^  = 2 TT  K p [(r2  + zi^)V2  - zi]  (2-32) 

Then,  for  Zj  = 0,  Zj  = »,  equation  2-29  further  reduces  to: 

= 2 TT  K p R (2-33) 

2.1.10  Vertical  Right  Circular  Cylinder,  General  Expression 

The  coordinate  system  \ised  to  define  the  cylinder  is  shown  in 
Figure  2-16.  The  gravitational  effect  of  the  cylinder  is  then  given  by 
equation  2-3**. 

g^  = K p R Gj  (2-31*) 

The  function  G5  is  a very  complex  hyperbolic  sine  function  which  is 
impractical  to  evaluate. 

Figure  2-17  gives  Gs  in  template  form  as  a function  of  the 


ratios  x/R  and  R/z.  The  template  is  taken  from  Talwani  (1973),  and  is 
reproduced  with  the  permission  of  the  author  and  publisher. 


The  function  G5  is  given  for  an  outcropping  vertical  circvilar 
cylinder  of  radius  R,  depth  zz  with  a vertical  axis  at  distance  x away 
from  the  computation  point.  The  effect  of  a cylinder  lying  between 
depth  z\  and  zz  is  then  obtained  by  subtraction. 


®Vi  2 ®V2  “ ®Vi 


(2-35) 


where : 


g = The  gravitational  effect  of  a cylinder  between  z,  and  z,. 

^1,2  * ^ 

g^  = The  gravitational  effect  of  an  outcropping  cylinder  of  depth 


g^  = The  gravitational  effect  of  an  outcropping  cylinder  of  depth 


2.1.11  Right  Circular  Cone.  Point  on  the  Axis 


Figiire  2-18  shows  the  right  circular  cone  and  the  coordinate 
system  in  which  it  is  defined.  Equation  2-36  then  gives  the  gravitational 
effect  of  a right  circular  cone. 


g^  = 2 IT  K P ^Z2  - Z^  - cos ^3  [(Z2^  + R2^)*/^  - Z^]  - 

(Z2^  + R2^)^'^^  + (z2  “ Z ) sec3  + z cos3 

*0 ^ nrrssisf ^ — 

o 


(2-36) 


If  z =0,  equation  2-36  is  simplified  as  follows: 
o 


gular  coordinate  system  as  shown  in  Figure  2-19*  The  gravitational  effect 
of  cylindrical  sectors  and  compairtments  Eire  given  by  equations  2-Uo  and 
2-kl  respectively. 


g^  = K a p [h  - (r^  - rs)] 

g^  = K a p (rj  + ra  - ri  - r4) 

where : 

di  = (xi^  + yi^)'/^ 

da  = (xa^  + 

n = (di^  + 


Ik. 


U2 


(Sector)  (2-Uo) 

(Compartment)  (2-Ul) 

(2-U2) 


Figure  2-19 


Cylindrical  Sectors  and 
Compartments 


rz  = rs 


= {dz^  + 


rs  = [di^  + (h  + z)2]V2 
= [da^  + (h  + z)^]^/^ 


(2-1^2) 

cont. 


Right  circular  cylinders  are  used  to  model  igneous  plugs  and 
intrusions  of  finite  cross-sectional  areas,  whereas,  right  circular  cones 
can  be  used  to  model  topographic  features.  Conical  sections,  and  cylin- 
drical sectors  and  compartments  form  the  basis  of  some  computation  schemes 
for  terrain  corrections . 

2.2  Computation  Schemes  for  Three-Dimensional  Bodies  of  Arbitrary  Shape 
The  computation  schemes  discussed  in  this  section  are  based  on 
equations  1-7,  1-8,  1-10,  and  1-12  in  the  section  on  mathematical  devel- 
opment. These  equations  define  the  triple  integration  in  various  coordinate 
systems,  as  restated  in  equations  2-43  to  2-46: 

Cartesian  Coordinates  (x,y,z): 


(2-43) 


Spherical  Coordinates  (r,4>,a): 


= K p 


sintfi  cos4)  dr  d<p  da 


(2-44) 


Cylindrical  Coordinates  (R,z,a): 


R z dz  dR  da 


^ ' j J J 


(2-45) 


Conical  Coordinates  (R,c|),a): 


Sv  = K P 


sincf)  dR  d(})  da 


(2-46) 


44 


where  the  parameters  x,  y,  z,  r,  R,  ()),  and  a are  defined  in  Figure  1-3. 
The  triple  integration  can  be  carried  out  in  several  ways. 

The  most  obvious  computation  schemes  involve  summations  over  volume 
elements,  horizontal  line  elements^or  vertical  line  elements.  The 
summation  over  the  volume  elements  is  a triple  niimerical  integration 
based  on  equation  2-2,  whereas,  the  summation  over  the  horizontal  or 
vertical  line  elements  involves  a single  analytical  integration  based 
on  either  equation  2-8  or  equation  2-11  and  a double  numerical 
integration.  Though  such  computation  schemes  are  arithmetically 
simple,  they  suffer  from  two  disadvantages.  First,  the  schemes  require 
the  body  be  divided  into  a large  number  of  volume  or  line  elements. 
Second,  the  effect  of  a unit  voliome  or  line  element  is  inversely  pro- 
portional to  its  distance  from  the  origin.  Therefore,  the  volume 
or  line  elements  must  be  of  variable  size  in  different  parts  of  the  body 
to  maintain  a uniform  degree  of  accuracy.  The  numerical  integrations 
in  volimie  or  line  element  computation  schemes  may  be  improved  by  replac- 
ing the  simple  summation  with  some  numerical  quadrature  technique  such 
as  Simpson's  Rule. 

The  most  commonly  used  computation  schemes  for  three-dimensional 
bodies  involve  various  combinations  of  analytical  and  numerical  inte- 
grations to  compute  the  combined  effects  of  horizontal  laminae,  vertical 
laminae,  cylindrical  wedges,  or  vertical  prisms. 

2.2.1  Schemes  Based  on  Rectangular  Vertical  Laminae 

The  three-dimensional  body  is  divided  into  a suitable 

1*5 


number  of  vertical  laminae.  Then,  each  lamina  in  the  plane  x = 
is  approximated  by  a suitable  polygon  ABC... A with  rectangular  corners 
and  thickness  Ax,  Figure  2-20.  The  gravitational  effect  of  lamina 
ABC... A is  then  stated  in  terms  of  G(y/x,  z/x)  , as  given  in  equation 
2-1+T. 


= K p Ax  [G(y./Xj,  a./Xj)  - GCy^^^/Xj,  . 


' 


(2-i+T) 


Thus , the  analytical  surface  integration  is  carried  out  in  y and  z 
with  the  numerical  integration  carried  out  along  x.  The  vertical 
lamina  method  is  applicable  to  horizontally  elongated  bodies,  but 
it  should  not  be  used  to  bodies  outcropping  near  the  computation  point. 

2.2.2  Schemes  Based  on  Rectangular  Horizontal  Laminae 
' The  three-dimensional  body  is  divided  into  a suitable 

number  of  horizontal  laminae.  Then,  an  arbitrary  lamina  in  the  z = Zj 
plane  is  approximated  by  the  polygon  ABC... A with  rectangular  corners 
and  thickness  Az,  Figure  2-21,  The  gravitational  effect  of  lamina  ABC... A 
is  then  stated  in  terms  of  G2(xy/xyz)  as  given  in  equation  2-U8. 


g^  = K p Ax  10i(x.y./x.y.Xj)  - Ox  ) 


* 02(x.^2yj^2/x.,2yi,2Xj) 1 

The  analytical  surface  integration  is  carried  out  in  x and  y and  the 
numerical  integration  carried  out  over  depth  z. 

2.2.3  Schemes  Based  on  Arbitrary  Horizontal  Polygonal  Laminae 
The  equations  presented  in  this  section  are  based  on  the 
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Figure  2-21 

Horizontal  Lamina,  XY  Plane 


Figure  2-22 


Arbitrary  Horizontal 
Polygonal  Lamina 
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original  derivations  by  Talwani  and  Ewing  (i960).  The  complete 

step-by-step  derivation  is  also  given  by  Clermont  (1967). 

The  body  is  approximated  by  a suitable  number  of  polygonal 

laminae.  A given  lamina  in  the  z = z.  plane  is  approximated  by  the 

J 

polygon  ABC... A of  thickness  Az,  Figure  2-22. 

The  gravitational  effect  of  the  lamina  ABC... A is  first 
expressed  as  a surface  integral  in  cylindrical  coordinates.  The  surface 
integral  is  then  transformed  into  a line  integral  over  the  polygonal 
lamina  boundary. 

The  gravitational  effect  of  the  triangular  lamina  P'BC  ‘ 
in  Figure  2-22  is  given  by  equation  2-hp. 


g^  = K p Az 


“i+1 


^ I 

J 


(2-1+9) 


The  parameter  r is  given  by  equation  2-50: 


(2-50) 


The  integral  in  equation  2-U9  is  solved  to  give  the  gravitational  effect 
of  the  triangular  lamina  P^BC.  The  gravitational  effect  of  the  entire 
polygonal  lamina  ABC... A is  obtained  by  a s;mimation  over  all  the  sides 
of  the  n-sided  polygon  as  given  in  equation  2-51. 


g = K p Az  [ 

i=l  L 


“i+1  - “i 


z.  cos3. 

(p"r+  z.^yi/^ 

^ J 


z cosy. 

* (p.“  t 8,”)'/' 


(2-51) 


Equation  2-51  is  rewritten  in  terms  of  x^,  y^,  and  ^i+i 

convenience  in  computer  application,  as  shown  in  equation  2-52. 


n 

g = K p ZSa  I 
i+1 


z q^  s 

s cos-  1 t . - sin-  1 i)T7-2 


z f.  s 
^ J 


(2-52) 


where : 


r.  = (x.2  + y.2)i/2 
1 1 1 

■•i.i.i  = “-i  - ^ ‘='i  ■ ^1*1’“'''“ 


, , , *i  *1.1  * >^1.1 

h = <“1.1  - “i> 


, ‘^1.1  *i  - ^1  ^.i>" 

Pi r~~= 


q.  = cosB.  = 
1 1 


i ,i+l 

i - vi^ " yj  ^yj  - yj^i^ 


r T 

i i,i+l 


(2-53) 


f.  = 
1 


Vi  ~ Vi^  * yj^i  ^^i  ~ yj^i^ 

^i+l  ^i,  i+1 


s = +1  if  sin  (0*^+-]^  ~ °^i^  ^ ^ 


s = -1  if  sin  (ot^+2^  ~ °*i^  ^ ^ 


The  total  gravitational  effect  is  again  obtained  by  niimerical 
integration  along  the  z eixis.  This  numerical  integration  can  be  performed 


the  rectangular  cross-section  ABCD  and  subtends  the  angle  ^ at  the  Z-axis 
as  shown  in  cylindrical  coordinates  in  Figure  2-23- 

The  gravitational  effects  of  cylindrical  sectors  and  compart- 
ments are  given  by  equations  2-kO  and  2-4l,  respectively.  Equations 
2-hO  and  2-hl  are  restated  in  cylindrical  coordinates  in  equations  2-5^ 
and  2-55. 

g^  = Aa  K p [(R2^  + - Zi  - (R2^  + + Z2)]  (Sector) 

(2-5U) 

and, 

g^  = Aa  K p [(R2"  + Zi2)1/2  . + 222)1/2  - (Ri2  + ^^2)1/2 

+ (Ri^  + Z2^)^/^]  (Compartment)  (2-55) 

The  numerical  integration  through  the  total  angle  a arovind  the  z-axis 
is  best  performed  by  some  numerical  technique  of  quadrature  such  as 
Simpson's  Rule. 
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Figure  2-25 

Cross  — Section  of  Arbitrary 
Cylindrical  Wedge 


the  RZ  plane  in  the  n-sided  polygon  as  shown  in  Figure  2-25- 

Talwani  (1973)  sixggests  such  wedge-shaped  compartments  be 
modeled  using  a variation  of  equation  2-39  which  gives  the  gravitational 
effect  of  a conical  section.  The  analytical  integrations  are  performed 
by  applying  equation  2-39  over  all  n-sides  of  the  polygon  A'S'C^D"*, 
summing,  and  then  multiplying  by  Aa  in  place  of  2ir.  The  result  is  an 
equation  of  the  form  of  equation  2-56. 


= AaK  p f^^i+1^  ^i+1^^  ■ ^^i^ 


+ (R.  cosB.  ^ sinB-  - z.  cosB-  sin^B- ) 


+ R^^)  cosB^  + + R^  sinB^  cosB^  - sin^B^  | 


(2-56) 


;rtiere: 


^^“^i  = ~~T. 


(^1  - (Vi  - 


zn.  = z.  - m.  R.  = Z..T  - m.  R. 

“i  1 11  1+1  1 1+1 


(2-57) 


r 


Ji  ,. 


by  Takin  and  Taiwan!  (1966). 

2.2-5  Schemes  Based  on  Right  Rectangulaj  Prisms 

The  three-dimensional  body  is  divided  into  a suitable  number 

of  rectangular  parallelepipeds  or  prisms.  The  gravitational  effect 

of  each  prism  is  computed  by  an  exact  expression  such  as  equation  2-58. 

The  total  gravitational  effect  of  the  three-dimensional  body  is  then 

the  algebraic  sum  of  the  prisms  making  up  the  body. 

The  prism  Po,  in  Figure  2-26  is  defined  by  opposite  corners 

P (0,0 ,0)  and  o (x,y,z).  The  gravitational  effect  of  Po  is  given  by 
1 1 1 

equation  2-58,  expressed  in  rectangular  coordinates. 


g^  = K p 


zi 


yi 


( 

\^2  ■ 


- sin 


_i 


zi 


xi 


(yi^  + Tirr^~*~y7T^ 


A [(xi^  + yi^  + - yj  [(xi^  + yi^)^/^  + yi] 

[(xi^  + yx'  + y\]  [(x,^  + . yj 


yi 


[(xi^  + yi^  + zi^)^/^  - xj  [(xi^  + yi^)^/^  + xj 
' 2 ^ [(xi^  + yi^  + + xj]  [ (xi^  + - xi]| 


(2-58) 


The  complex  functions  in  brackets  are  rewritten  in  functional  notation 
as  follows : 


1 

L 


57 


Then  equation  2-58  is  simplified  to  the  form  of  equation  2-60: 

= K p zi  G3(yi/xi,  zi/xi)  -(2- 

where  the  function  G3  is  evaluated  directly  from  Figure  2-27  for  given 
ratios  of  x,  y and  z. 

The  gravitational  effect  of  an  arbitrary  rectangular  prism 
oh  is  computed  by  adding  and  subtracting  the  gravitational  effects  of 
eight  prisms  as  shown  by  equation  2-6l. 


®oh  ■ ®Ph  ” ®Pf  * ®Pe  ” ®PJ  " ®Pc  ®Pb  ■ ®Po  * ®Pd 
The  subscripts  denote  opposite  corners  of  the  prisms  in  Figure  2-26. 

Nagy  (1966)  presented  refined  versions  of  equation  2-58 
suitable  for  modeling  and  terrain  effect  computations.  These  equations 
are  advantageous  because  they  contain  none  of  the  inherent  errors  asso- 
ciated with  approximate  formulas.  Also,  the  Nagy  equations  are  readily 
adaptable  to  computer  applications. 


3.  TWO-DIMENSIONAL  ATTRACTION  FORMULAS 

3.1  Some  Simple  Two-Dimensional  Attraction  Formulas 

The  equations  presented  in  this  section  are  derived  from  equation 

1-llt  in  the  section  on  mathematicsil  development.  The  two-dimensional 

bodies  are  defined  in  a rectangular  coordinate  system  in  x and  z. 

The  angles  d)  are  meaSTired  clockwise  from  the  x-axis. 
n 

3.1.1  Infinite  Horizontal  Rectangular  Prism 

Figure  3-1  gives  the  geometrical  and  angular  relationships 
of  an  infinite  horizontal  rectangular  prism.  The  gravitational  effect 
is  then  given  by  equation  3-1. 


g = 2 K p ■X2  In  — - Xj  In  — + Z2  {())2  " (fru)  “ Zi  ((|)i  ” 4)3) 
* , ^3  -^1 


(3-1) 


where : 
^1 


= (xi^  + 

= {X22  + zi^)V2 
= (X2^  + Z2^)^/^ 


(3-2) 


<t>l 


_ TT  * -1 

- - - tan  ' — 
2 zi 


TT 

4)2  = 2 - tan-1  _ 


^2 

A - 4.-1 

‘*’3  = ? - IT 


X2 

4.4  = ^ - tan- 
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Figure  3-2 


Infinite  Horizontal  Ri 
Rectangular  Pris 


The  infinite  right  triangular  prism  is  used  to  model  bodies  which  are 
genereilly  asymmetrical  in  cross-section. 


3.1.3  Symmetrical  Anticline  of  Infinite  Extent 

The  symmetrical  anticline  is  formed  by  the  sum  of  two  right 
triangular  prisms  of  infinite  extent  as  shown  in  Figure  3-3.  The 
gravitational  effect  is  given  by  eq.uation  3-6. 


g^  = 2 K p [xj  sin  i] 


sin  i In  — + cos  i ( ())2  + ({>3  “ 2(|)j) 


1-2 

sin  i In  ~ * cos  i (4)2  - 4)3) 


+ Z2  (412  “ 4)3)/  (3-6) 


where : 
r 


r,  = 


r,  = 


(Xj2  + Z j2)  1/2 

[(Xj  - a)2  + Z22]i/2 
[(x,  + a)2  + Z22]i/2 


= JI_  tan-1  — 


•fi  = 


4)2  = f - tan"i 


(3-T) 


4,3  = ^ - tan 


- 1 


Xj  + a 


cos  i = ip-r^ynT 
sin  i = 
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Geldaxt,  Gill,  and  Sharma  (1966)  give  an  excellent  discussion  of  gravity 
anomalies  of  two-dimensional  faults. 

3.1.6  Inclined  Prism  of  Infinite  Extent  ^ 

The  inclined  bed  of  infinite  extent  is  formed  by  the  dif- 
ference of  two  offset  slopes  as  shown  in  Figure  3-6.  Equation  3-ll* 


gives  the  gravitational  effect. 

[ ra  ra 

= 2 K p [xi  sin  i + zi  cos  i]  [sin  i In 

+ cos  i (4)2  - 4>i  + <{>3  - <}'4)]  + a sin  i [in  — + cos  i (4)4 

+ Z2  (4)2  - 4)4)  - zi  (4)1  - <))3)|^ 
where : 


The  infinite  inclined  prism  is  used  to  model  inclined  dikes  or  similar 
structures. 


3.1.7  Two-Dimensional  Criterion 

There  are  no  geologic  or  geophysical  structirres  that  are 
truly  two-dimensional.  However,  structures  do  exist  that  are  sufficiently 
elongated  in  a single  horizontal  direction  to  be  considered  two-dimensional. 
The  use  of  two-dimensional  approximations  in  modeling  and  interpretation 
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Figure  3-6 


Inclined  Prism  of  Infinite  Extent 
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is  based  on  the  fact  that  the  gravitational  effect  of  more  distant 
masses  becomes  negligible.  It  is  important  to  estimate  the  error 
resulting  from  making  the  assumption  of  two-dimensionality. 

The  maximum  error  resulting  from  a two-dimensional  assump- 
tion can  be  estimated  by  substituting  a finite  horizontal  line  element 
for  the  two-dimensional  body  as  shown  in  Figure  3-7-  The  distance  r 
is  measured,  in  the  xz  plane,  from  the  computation  point  or  origin  to 
the  point  on  the  cross-section  of  the  two-dimensional  body  farthest  from 
the  origin. 

r = [ (x  + Ax)^  + (z  + Az)^]V^  (3-16) 

If  the  cross-sectional  area  of  the  body  is  small.  Ax  and  Ay  will  be 
negligible.  Then, .r  becomes  the  distance  in  the  xz  plane  from  the 
origin  to  the  line  element  and  is  computed  from  equation  3-17. 


r = (x^  + 


(3-17) 


Then,  for  a given  point  on  the  X-aocis , r is  dependent  on  z and  the 
maximum  percentage  error  p,  in  making  the  two-dimensional  assumption 
is  then  estimated  by  equation  3-l8: 

p = (1  - e)  100^  (3-18) 


where : 

^ " (r^  ; ^2)172-  (3-19) 

Table  3-1  gives  the  values  of  p for  various  ratios  of  y and  r. 
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TABLE  3-1 


Two-Dimensional  Criterion 


y 0 .5r  r 1.5r'  2r  3r  Ur'  5r  lOr  “> 

pS?  100  55.3  29.3  16.8  10.6  5.1  3.0  1.9  0.5  0 


3. 2 Computation  Schemes  for  Two-Dimensional  Bodies  of  Arbitrary  Shape 
The  equations  presented  in  this  section  form  the  basis  for  many 
automated  modeling  and  inversion  schemes.  Such  schemes  are  used  to 
model  two-dimensional  bodies  of  arbitrary  shapes  and  variable  density 
contrast. 

3.2.1  Equations  for  Surface  Integration' 

The  gravitational  effect  of  a two-dimensional  body  of  arbi- 
trary cross-sectional  area  A is  obtained  by  integrating  the  effects 

s 

of  infinite  horizontal  line  elements  of  area  dA  over  A , Figure  3-8. 

s 

The  general  form  of  the  integral  is  restated  in  equation  3-20. 


= 2 K p 


(3-20) 


Equation  3-20  is  expressed  in  different  coordinate  systems  as  follows: 
Cartesian  Coordinates,  (x,  z): 


- „ f f z dx  dz 

ev  = 2 K P J J (P-V  --rj 


(3-21) 


Polar  Coordinates,  (r,0): 


g •=  2 K p sin0  d0  dr 

s 


(3-22) 


Angular  Coordinates,  (x,6); 


g = 2 K p 1 tan0  d6  dx 
V /a  i 


(3-23) 


Angular  Coordinates,  (z,6); 


g^  = 2 K p J j de  dz 
^s 


(3-2lt) 


The  above  equations  can  be  evaluated  by  a simple  summation  of  the  effects 

of  individual  horizontal  line  elements  over  the  area  A . However,  the 

s 

gravitational  effect  of  a line  element  increases  very  rapidly  as  its 
disteince  from  the  origin  decreases. 

3.2.2  Equations  in  Terms  of  Line  Integrals 

It  is  more  practical  and  convenient  to  solve  the  two-dimen- 
sional problem  in  terms  of  line  integrals  rather  than  surface  integrals, 
Hubbert  (I9lt8).  Equations  3-22,  3-23  and  3-24  are  rewritten  as  line 
integrals  in  equations  3-25  through  3-28  where  the  line  integral  L 


is  in  the  cross-sectional  plane  of  the  body. 

For  the  polar  coordinates,  (r,6)  figure  3-9: 


g = 2 K p Ar  sin0  d0 

=>v  j 


= 2 K p Ar  [(cos02  - cos0i)  + (cos04  - cos03)  + ...] 
s 


For  the  angular  coordinates,  (x,0)  figure  3-10: 


g^  = 2 K p Ax  I tan9  d0 


= 2 K p Ax  In 


cos  02 

COS  01 


+ In 


COS0^ 

COS03 


n 

= 2 K p Ax  I 
i=l 


cos0.^ll 


COS  8. 
1 


(3-26) 


Equation  3-2k  is  the  simplest  and  most  convenient  to  xise. 


The  integral  can  be  expressed  in  terms  of  Az  or  A0  as  is  shown  in  equa- 
tions 3-27  and  3-28. 

For  the  angular  coordinates,  (z,0)  figures  (3-ll)  and  (3-12): 

= 2 K p Az  I d0 
L 


= 2 K p Az  [(02  - 0i)  + (04  - 63)  + ...] 
n 

= 2 K p Az  [ (0.^^  - 0.)  (3-27) 

i=l 

Or,  by  changing  the  order  of  integration,  equation  3-2^+  becomes: 

g^  = 2 K p A0  I dz 
L 

= 2 K p A0  [(z2  - Zi)  + (z4  - Z3)  + ...] 


n 

= 2 K p A0  [ (z.^^  - z. ) (3-28) 

i=l 

Any  of  the  integration  schemes  defined  by  equations  3-25 
through  3-28  can  be  used  to  numerically  evaluate  the  line  integral  L. 
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However,  the  most  commonly  used  schemes  are  based  on  the  approach  of 
Taiwan!,  Worzel  and  Landisman,  (.1959)  • It  is  a practical  and  convenient 
approach  in  which  the  two-dimensional  cross-section  is  approximated  by 
an  n-sided  polygon. 

3.2.3  Gravitational  Effect  of  an  n-Sided  Polygon 

The  gravitational  effect  of  an  arbitrary  two-dimensional  body 
of  triangular  cross-section  PCD,  Figure  3-13,  is  given  by  equation  3-29. 


= 2 K p 


r®i.i 

0. 

1 


z d6 


(3-29) 


The  gravitational  effect  of  the  polygon  ABCD...A  is  obtained  by  inte- 
grating equation  3-29  for  each  triangular  cross-section  and  summing  the 
results,  as  in  equation  3-30. 


N 1 X.  z...  - z.  x.,^ 

_ „ ^ „ r 1 1+1  1 1+1 

“ 2 K p I -7-  p + ( z P" 

i=l  r^i  ''i+l^  ^^i  ^i+1^ 


l<Vl  - - Vl>  <»1.1  - »l”('  <3-30) 

1 


B3 


:'i 


k.  DENSITY  DETERMINATION 
U.l  Density  and  Porosity  Defined 

Density  is  defined  as  the  ratio  of  mass  to  volume.  If  a given  rock 
sample  is  100^  solid  state  rock  material,  density  is  defined  by  equation 


(W) 


In  reality,  rock  specimens  are  composed  of  solid  state  rock  mauerial, 
and  liquid,  and/or  air  filled  pore  space.  Therefore,  equation  U-1  must 


be  rewritten  as : 


(Ml  + M2  + M3) 
P " (Vi  + V2  +”V3) 


(U-2) 


where : 


Mj , M2,  M3  = The  masses  of  the  solid  state,  liquid,  and  air  spaces, 
respectively. 

^2  > ^^3  “ volumes  of  the  solid  state,  liquid,  and  air  spaces, 
respectively. 

Equation  1-2  points  out  the  inherent  ambiguity  in  density  determination. 
Therefore,  in  practice  it  is  necessary  to  define  density  in  terms  of 


some  limiting  relationships. 

Precise  laboratory  measurements  yield  accurate  values  of  the  dry 
and  saturated  densities  and  as  defined  in  equations  1-3  and  1-1: 


Ml  + M3 
•^d  " Vi  + V3 
Ml  + M2 
^s  ""  Vi  + Vz 


(1-3) 


(1-1) 
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Porosity  p is  the  ratio  of  the  volume  of  liquid  and  air  space  in  a given  \ 

i 

rock  sample  to  the  total  volume  as  given  in  equation  4-5.  i 

V2  + V3 

" Vj  + Vj  + Vj 

The  true  rock  sample  density  depends  on  mineralogic  composition, 
porosity  and  degree  of  saturation  and  is  defined  by  the  following 
range : 

Pd  < Pv  1 Ps 

1 

Ideally,  large  numbers  of  rock  samples  are  measiored  and  processed  statis-  ^ 

tically  to  yield  a best  estimate  of  the  tri’^  density  of  the  rock  samples.  j 

i 

4. 2 Generailized  Density  Kelp.tionships 

Generalized  density  rr lationships  are  best  discussed  in  terms 
of  rock  types  (i.e.  , igneous,  metamorphic  and  sedimentary).  The  density 

of  igneous  and  metamorphic  rocks  depends  mainly  on  mineralogic  composi-  i 

i 

tion  because  the  porosity  of  these  rocks  is  usually  less  than  2!?  - k%.  i 

The  density  relationships  in  sedimentary  rocks  are  highly  variable. 

Carbonate  rocks  containing  few  or  no  solution  cavities  are  the  only 
sedimenteury  rocks  in  which  mineralogic  composition  is  generally  more 
important  than  porosity.  Pick,  Picha  and  Vyskocil  (1973)  estimate  the 
following  porosity  ranges  of  sedimentary  rocks,  based  on  the  degree  of 

J 

consolidation.  Unconsolidated  sedimentary  rocks  have  porosities  ranging 
from  2 55^  to  90%»  Sedimentary  rocks  consolidated  by  diagenesis  (i.e., 
static  processes  of  lithification)  have  porosities  tending  towards  185?, 
whereas,  sedimentary  rocks  consolidated  by  severe  orogenesis  (i.e., 
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dynamic  processes  of  lithification)  have  porosities  averaging  about  k%. 

In  general,  the  density  of  sedimentary  rocks  increases  with  depth  of 
burial.  However,  this  general  relationship  is  variable  and  must  be 
determined  locally. 

Residual  gravity  anomalies  are  caused  by  lateral  density  or  mass 
variations  in  the  earth's  crust  and  upper  mantle.  If  the  earth  were 
composed  of  material  in  layers  of  laterally  uniform  mass,  there  would 
be  no  residual  gravity  anomalies  no  matter  what  the  vertical  veuriation. 

Many  gravimetric  investigations  are  primarily  dependent  on  a knowledge 
of  mass  distribution  and/or  local  isostatic  conditions.  However,  mass 
distributions  may  be  diffj.cult  to  determine  and  express  due  to  some 
inherent  ambiguities.  Regional-residual  sepfration  techniques  may  not 
realistically  separate  the  regional  and  residual  fields.  Also,  mass 
distributions  may  be  subtJ.e  continuous  functions  of  position  rather  than 
clearly  defined  discrete  fimctions  of  position.  Such  ambiguities  can 
significantly  affect  the  results  of  gravitational  modeling  and  interpre- 
tation. 

An  examination  of  the  general  integral  equations  for  the  gravita- 
tional effect  of  two  or  three-dimensional  bodies  shows  that  density  is 
independent  of  the  geometry  of  the  body.  In  fact,  the  other  unknown 
parameters,  depth  and  size  are  indirectly  dependent  on  the  initial  density 
assumptions.  Therefore,  the  initial  density  approximations  must  be 
determined  by  methods  independent  of  geometry  such  as  the  laboratory 
methods  already  mentioned.  However,  densities  determined  in  the  laboratory 
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may  not  be  representative  of  large  masses  of  rocks  in  situ, 

U. 3 Nettleton's  Method  of  Density  Profiling 

Nettleton  (1939)  describes  a method  of  determining  "effective" 
densities  from  gravity  measurements.  The  effective  density  is  the 
in  situ  density  of  leirge  topographic  masses  due  to  a variety  of  influ- 
ences such  as:  mineralogic  composition,  porosity  changes  in  rock  type, 
and  degree  of  isostatic  compensation. 

Gravity  and  elevation  meas\irements  are  made  along  a profile  usually 
perpendicular  to  the  topographic  structure-  Free-air  gravity  anomalies 
are  then  computed  at  n observation  points  j=l,  n,  using  equation  U-7:' 


.■"obs'  "j 

i J 


(M) 


where : 


= Free-air  gravity  anomaly  at  point 

'j 

j®obs)  " Observed  gravity  at  the  J point 

J 

th 

Yj  = Theoretical  gravity  at  the  point 

h.  = Elevation  at  the  point 

J 

iSg^  = Free-air  reduction  (0.3086  mgals/meter,  elevation  in  meters) 


A geophysically  realistic  range  of  densities  is  chosen,  pj  to  and 
a Bouguer  gravity  anomaly  profile  is  computed  for  each  density  using 


equation  U-8. 


(U-8) 


where : 

27r  K p h = The  Bouguer  reduction  for  the  density  at  the 
observation  point,  i=l,  m. 

If  the  elevation  range  is  leirge,  the  Free-air  gravity  anomalies  should 
be  terrain  corrected  for  the  surrounding  topography.  The  resiilting  Bouguer 
gravity  anomaly  profiles  are  plotted  against  distance  and  compared  with 
elevation,  as  in  Figure  it-1. 

Equation  k-8  shows  the  Boiaguer  gravity  anomalies  to  be  directly 
proportional  to  elevation  for  all  densities  less  than  the  true  density 
(i.e.  , < p^)  and  the  inverse  holds  for  (p^  > p^).  Thus,  the  Bouguer 

gravity  anomaly  profile  for  the  correct  density  is  the  one  reflecting 
the  least  dependency  on  elevation  (i.e.  , the  profile  should  be  nearly 
flat  directly  over  the  topographic  feature). 

Figxore  4-1  illustrates  the  method  as  applied  over  a hypothetical 
seamount.  Bouguer  gravity  anomaly  profiles  are  plotted  for  four  densities. 
The  profiles  for  p = 1.7  gm/cm^  and  p = 2.2  gm/cm®  clearly  show  a direct 
relationship  with  elevation,  whereas  the  profile  for  p = 2. 4 gm/cm® 
shows  an  inverse  relationship  with  elevation.  Thus,  the  effective  density 
is  Py  = 2.3  gm/cm®  which  agrees  with  the  density  used  in  constructing 
the  hypothetical  seamovint  model. 


's  Method  of  Density  Determination 


Jxuig  (1943)  expresses  the  Nettleton  method  in  mathematical  terms 
that  can  be  applied  to  surface  distributions  of  points  as  well  as  profiles. 
Gravity  euid  elevation  measiurements  are  made  over  the  area  under  study. 


I 


Equation  U-7  is  then  used  to  con5)ute  Free-air  gravity  anomalies  at  the 
n-observation  points  suid  the  Bouguer  gravity  anomalies  are  conq)uted  for 
one  of  the  densities  (i=l»  m)  in  the  range  of  geophysically  realistic 
densities. 

The  method  of  least  squares  regression  analysis  is  used  to  determine 
effective  density.  The  Bouguer  gravity  anomalies  for  each  of  the  given 
densities  can  be  expressed  as  a function  of  elevation  in  the  form  of  a 


regression  line  by  equation  U-9: 


•;^biij  = "i " ^ 


(U-9) 


where: 


’ ^g,  = The  Bouguer  gravity  anomaly  for  the  i^^  density  p.  at 

^ th 

the  J point 

til 

a^  = The  intercept  of  the  regression  line  for  the  i density  p^ 

b^  = The  slope  constant  of  the  regression  line  for  the  i^^ 

density  p^^ 

First,  the  correlation  coefficient,  r^,  is  computed  for  each  regression 
line  using  equation  U-IO; 


1=1  . L 


(U-10) 


Ah 

^ J=1  . ij 


where : 
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i [ 


1^1.  = 


(1+-11) 


n 

h = 

n 


If  the  correlation  coefficient,  r^  is  equal  to  zero,  the  anomalies  j 

are  uncorrelateci  for  a given  density  p^.  Thus,  the  density  is  the 
real  effective  density  of  the  topographic  mass. 

If  the  correlation  coefficient  is  not  equal  to  or  sufficiently  close 
to  zero,  the  density  p^  is  not  the  effective  density.  Then,  for  r^  f o, 
the  slope  constant  of  the  regression  line  is  computed  by  equation  U-12. 


(i*-12) 


I 


Equation  4-9  is  then  rewritten  in  the  following  form: 


(4-13) 


where; 


= The  Bouguer  gravity  anomaly  computed  using  the  true  effec- 
tive  density  at  the  computation  point 


” The  true  effective  density. 
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The  slope  constant  in  equation  U-13  is  defined  by  the  following 


equation: 


b.  = 27r  K (p.  - P.  ) 
1 t 1 


(1|-1U) 


Then,  the  correct  value  of  the  effective  density  P is  con5)uted  by 


equation  U-15: 


b. 

^t  ^ ^i  2inc 


('+-15] 


where,  b^  has  been  computed  by  equation  4-12. 

The  Nettleton  aoid  Jung  methods  both  give  the  effective  densities 
of  the  topographic  structures  or  masses  between  the  highest  and  lowest 
elevations  over  the  profile  or  area.  The  methods  do  not  work  well  in 
areas  of  low  relief,  areas  of  complex  lateral  density  variations  or  areas 
of  considerable  basement  relief  underlying  the  sedimenteiry  layers. 


5.  SOME  SIMPLE  TECHNIQUES  FOR  DEPTH  DETERMINATION 
The  purpose  of  this  section  is  to  present  some  simple  techniques  and 
formulas  for  the  determination  of  depth  and  possible  shapes  and  sizes. 
There  is  no  xmique  mathematical  solution  to  the  determination  of  the 
mass  distribution  of  the  source  of  a residual  gravity  anomaly  field. 
Therefore,  some  realistic  assumptions  and/or  estimations  must  be  made 
regarding  the  unknown  parameters:  depth,  size  and  shape.  Such  ass\anp- 
tions  are  generalizations  but  yield  realistic  first  approximations  when 
detailed  knowledge  of  structure  is  lacking. 

Cone  of  Solutions  Defined 

Nettleton  (1971)  suggests  the  *'Cone  of  Solutions"  as  an  initial 
approach  to  the  inverse  problem.  The  "Cone  of  Solutions"  states  that 
the  spherical  or  point  mass  of  given  density  is  the  deepest  possible 
singular  mass  configuration  that  can  cause  a given  residual  gravity 
anomaly  field.  However,  other  lenticular  or  rectangular  mass  configtura- 
tions  with  that  density  are  possible  at  shallower  depths,  as  shown  in 
Figiire  5-1.  All  possibilities  in  the  "Cone  of  Solutions"  should  preserve 
mass  (i.e.,  the  product  of  volume  and  density  contrast  must  be  the  same). 
The  following  paragraphs  describe  some  inversion  techniques  based  on 
some  possible  simple  mass  configurations  in  the  "Cone  of  Solutions." 

Each  technique  is  then  illustrated  using  the  hypothetical  residual  gravity 
anomaly  profile  shown  in  Figure  5-2. 

Skeels  (1963)  shows  that  the  deepest  geophysically  feasible  mass 
configuration  must  be  a "vertical-sided  mass"  (i.e.,  infinite  prism 
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Figure  5-2 
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or  right  circiilar  cylinder  for  two-  and  three-dimensional  cases , respec- 
tively). He  presents  empirical  charts  for  the  rapid  determination  of 
depths  and  widths  of  prisms  or  cylinders  for  a given  density  contrast. 
5.2  Spherical  Mass  Configuration 

The  gravitational  effect  of  the  spherical  mass , shown  in  cross- 
section  in  Figiire  5“3>  is  given  by  equation  5-1  = 


= I TT  K p Z 


(5-1) 


The  depth  to  the  center  of  mass,  z,  and  radius  of  the  sphere,  R,  are 

determined  by  the  procedirre  outlined  below. 

A profile  of  g^  versus  x is  constructed  through  the  maximum  residual 

gravity  anomaly,  as  shown  in  Figure  5~3.  The  points  A and  A'  on  the  x-axis 

correspond  to  g = 1/2  g^  , where  is  the  maximum  value  on  the  residual 

o o 

gravity  anomaly  profile;  and  the  distances  1oa|  = loA"!  = Xi/2.  Then, 
the  depth  to  the  center  of  mass  is  given  by  equation  5-2: 


z = 


Xl/2 

XPT^~Tp7^\ 


= 1.305  y^i/2 


(5-2) 


If  the  density  or  density  contrast  is  known,  the  radius  R is  given  by 
equation  5-3. 


9 « 2 

3z  g„ 


R = 


4ti  K p 


11/3 


(5-3) 


The  "eqxxivalent  sphere"  of  a given  density  or  density  contrast  is 
the  spherical  mass  at  depth  z which  will  most  nearly  approximate  the 
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residual  gravity  anomaly  profile.  In  the  model  illustrated  in  Figtire  5-3, 

the  density  is  assumed  to  be  0.5  gm/cm®.  The  "half -width" , Xi/2,  and 

maximum  residual  gravity  anomaly,  g , are  measured  from  the  profile: 

o 

xi/2  = 3.07  km  and  g =7  mgals.  Equations  5-2  and  5-3  yield  the  follow- 

'^o 

ing  values  for  z and  R:  z = U.OO  km,  R = 2.00  km. 

5.3  Mass  Distributed  Along  an  Infinite  Horizontal  Cylinder: 


The  gravitationed.  effect  of  the  infinite  horizontal  cylinder  shown 
in  cross-section  in  Figure  3-k  is  given  by  equation  5-^- 

_ _ 2 IT  K p z r1 


®v  + z* 


(5-U) 


The  profile  of  g versus  x is  constructed  such  that  |0A|  = |0A^1  = Xi/2« 


Then,  the  depth  to  the  center  of  mass  is  given  by  equation  5-5: 

z = Xi/2 


(5-5) 


The  density  or  density  contrast  is  then  used  to  compute  the  radius  by 


equation  5-6: 


2 IT  K p 


(5-6) 


With  the  same  values  for  p , Xiy2  and  g^  , equations  5-5  and  5-6  yield 

o 

the  following  values  of  z and  R for  the  equivalent  infinite  horizontal 
cylinder:  z = 3.07  km  and  R = 1.01  km. 

5.4  Mass  Distributed  Over  an  Infinite  Horizontal  Plane: 

The  infinite  horizontal  plane  may  be  substituted  for  a very  thin 
two-dimensional  rectangular  prism  as  shown  in  Figure  5-5.  The  approxima- 
tion is  valid  provided  the  thickness  of  the  prism  t is  small  compared 
with  its  width  w (i.e.  , t«w). 
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Infinite  Cylindrical  Mass 


Figure  5-5 





The  surface  density  constant,  y,  of  the  infinite  plane  is  defined  in 
terms  of  the  surface  density  of  the  two-dimensional  prism  by  equation  5-7. 


y = p (z2  - z i)  (5-T) 

The  gravitational  effect  of  the  infinite  plane  is  defined  by  equation  5-8: 

r 


. 2 K p itao-.  1^1  - tan-.  (i| 

= 2 K y ((,  (5-8) 

where  x is  the  distance  from  the  computation  point  to  the  near  edge  of 
the  infinite  plane. 

The  points  A,  A'’,  and  B and  on  the  x-axis  correspond  to  g^  = (1/^ 
g^  and  g^  = (l/l()  g^  , respectively.  Then,  the  distances  |0A|  = = 


o o 

Xj/j  and  |OBj  = |0B^|  = The  depth  to  the  infinite  plane  is  then 


given  by  equation  3~9- 


(Xi/4  - Xj/2) 


z = 


(5-9) 


2x1/2 

Then,  the  width,  W,  of  the  plane  is  determined  by  equation  5-10. 

W = 2 (x^/2  - z^)^/^  (5-10) 

Equation  5-7  can  be  used  to  compute  y only  if  values  are  known  for  Zi  and 
Z2.  However,  if  zi  and  zz  aren't  known,  y,  can  be  evaluated  from  equation 


5-8  using  z and  g as  shown  in  equation  5-10: 
o 


y = 


2 K ({)  max 


(5-11) 


1 
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With  X 2 s=  3. 07  km  and  = 7 mgals,  equations  5-9,  5-10  and  5-11 


o 

yield  the  following  values  for  z,  W and  y for  the  eqxii valent  horizontal 
infinite  plane:  z = 2.53  km,  W = 3. 1+49  km,  y=  .1*3973  gm/cm  ® (km). 

The  spherical  mass  amd  infinite  plane  are  the  two  limiting  cases 
in  the  "Cone  of  Solutions",  whereas  the  infinite  horizontal  cylinder 
is  an  intermediate  case  for  an  assumed  density  or  density  contrast. 

Thus  the  center  of  mass  of  the  mass  configuration  causing  the  residual 
gravity  anomaly  profile  in  Figure  5-2  must  be  at  a depth  between  2.53  km 
and  1*  km. 

The  methods  just  discussed  apply  to  symmetrical  residual  gravity 
anomaly  profiles  caused  by  either  spherical  or  two-dimensional  mass 
distributions.  Asymmetrical  residual  gravity  anomaly  profiles  must  be 
handled  by  more  refined  techniques  such  as  the  characteristic  curve 
method  described  by  Grant  and  West  (1965). 

5*  5 Mass  Distributed  Over  an  Infinite  Horizontal  Half-Plane: 

The  infinite  horizontal  half-plane  may  be  used  to  approximate  a 
vertical  fault  as  shown  in  Figure  5-6.  The  thickness  of  the  half-plane, 
is  then  the  vertical  displacement  along  the  fault.  The  one-limbed 
residual  gravity  euiomaly  profile  shown  in  Figure  5-7  is  characteristic 
of  a vertical  fault. 


Figure  5-  6 

ass  Over  an  Infinite 


6v  = 2 K y||-  tan-i^ 


(5-12) 


where  x is  the  distance  from  the  computation  point  to  the  ed^e  of  the 
infinite  half -plane. 

The  point  0 on  the.x-axis  must  be  determined  graphically  because 

it  corresponds  to  the  inflection  point  on  the  profile  directly  above 

the  edge  of  the  half-plane.  The  points  A and  B on  the  x-axis  correspond 

to  g = (1/2)  g and  g,  = ( 3/2 ) g , respectively , and  the  depth  z 
^o  B o 

is  given  by  equation  5-13. 


|0A|  = |0B|  s z 


(5-13) 


Eq\xation  5-12  is  then  evaluated  at  the  inflection  point  to  compute  a 


value  for  y,  as  shown  in  equation  5-1^: 


(5-14) 


Equation  5-13  and  5-14  yield  the  following  values  of  z and  y for  the 
residual  gravity  anomaly  profile  in  Figure  5“?!  z = 1 km,  y = .14906  x 
10®  gm/cm^. 

5.6  Estimation  of  the  Depth  to  the  Upper  Boimdary  of  a Body  of  Arbitrary 
Shape: 

The  general  solutions  to  the  problem  along  with  complete  derivations 
are  given  by  Bott  and  Smith  (1958).  The  equations  they  derived  are 
valid  in  a general  sense  and  are  applied  to  the  hypothetical  residual 
gravity  anomaly  profile. 

The  residual  gravity  anomaly  profile  referenced  is  given  in  Figure  5-2. 
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with  corresponding  and  x^  on  the  x-axis , 


Two  points  g and  g 

^2 


eure  chosen  such  that  g > g - Then,  A is  defined  by  equation  5-15 

"^1  ^2 


(A  > 1): 


A = 


(5-15) 


"V2 


The  maximum  possible  depth  to  the  upper  boundary  of  the  body  is  given 
by  equation  5-l6: 

z < ^5-16) 

Equation  5-lT  gives  z for  ni  > 1 where  m is  given  by  equation  5-l8: 


w 


1)V2 

The  parameter  w is  an  arbitrary  interval  on  the  x-axis. 


(5-lT) 


2g. 


m = 


’v(x) 


^®v{x+w)  ^ ®v(x-w)^ 

If  the  maximum  values  of  the  residual  gravity  anomaly  and  horizontal 
gradient  are  known,  the  depth  z is  given  by  equation  5-19: 


(5-18) 


z < .86 


^v^maxj_ 


ae„ 


(5-19) 


max 


The  same  derivations  yield  the  corresponding  equations  for  two- 
dimensional  bodies , with  m and  A already  defined. 


- A-1 


z < 


w 


(5-20) 

(5-21) 
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Z ^ .65 


'v(max) 


(5-22) 


max 


Equations  5-15  and  5-I6  yield  a value  of  z ^ 5-1  km  for  the  residual 
gravity  anomaly  profile  in  Figure  5-2»  where  |xj  - Xj  | = 5 km. 

5.7  Depth  to  a Density  Interface  Along  a Profile  With  Gravity  and 


Borehole  Information: 

The  density  interface  and  associated  residual  gravity  anomaly 
are  shown  in  Figure  5-8.  There  are  boreholes  along  the  profile  at  Xj* 
X3  and  X , and  the  interface  outcrops  at  Xi,such  that  the  depths  of  the 


interface  z(xi),  z(x2)  ...  z(x^)  are  known. 


The  interface  is  then  expressed  in  the  form  of  a series  as  in 
equation  5-23. 

z(x)  = ^2(gy)^  + . . . + A^(g^)’^ 

The  equations  are  set  up  as  in  equation  5 -2k. 

zi  = + Ai(g^^)  + A2(gy^)^  + . . . + 

Z2  = A^  + Ai{g  ) + A2(g  )"  + . . . + A^(g^  )"■ 


(5-23) 


(5-2k) 


^n  = Ao  + Ai{g^  ) + A2(g^  )^  + . . . + A^(g^  . 


The  solution  of  the  set  of  equations  ;yields  the  coefficients  Ai , A2  •••  A^. 
The  solution  is  unique  when  the  number  of  coefficients  equals  the  number 
of  observation  points  (i.e.  , i = n)  and  can  be  solved  as  simulteuieous 
equations. 
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However,  the  solution  is  not  tmique  when  the  nimber  of  observation  points 
is  greater  than  the  number  of  coefficients  (i.e.  , n > i).  The  best 
fit  solution  is  then  obtained  by  the  method  of  least  squares. 

If  the  depth  is  known  at  many  points  along  the  profile,  it  may 
suffice  to  graphically  interpolate  the  shape  between  the  known  depths , 
using  the  residual  gravity  anomaly  profile  as  control.  Both  the  graphical 
and  analytical  methods  are  independent  of  density. 


Anomalies  Alone: 

The  residual  gravity  anome.ly  profile  and  density  interface  are  shown 
in  Figure  5-9*  The  method  requires  that  the  densities  pj  and  pj  be 
known  or  reasonably  approxi'nsted.  The  density  contrast  is  defined  by 
equation  5~25: 

Pc  = Pi  - P2  (5-25) 


For  basin  structures,  p2  > pi  and  p^  is  negative. 

The  first  approximation  to  the  shape  of  the  interface  is  determined 
by  computing  depths  at  n-observation  points  along  the  x-axis  by 
equation  5-26; 


z. 

1 


2 IT  K p 

c 


(5-26) 


The  theoretical  gravity  anomalies  I g \ are  computed  for  the  current 

I \li 

depths  using  appropriate  equations,  i.e.,  Talwani's  equation  for  the 
gravitational  effect  of  two-dimensional  n-sided  polygon,  equation  3-31. 
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BL. 
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Gravity 


g 


The  theoretical  and  observed  residual  gravity  anomalies  are  compared 
and  differences  computed  using  equation  5-2T • 


= Sv.  - \]. 

1 k 1 


(5-2T) 


where : 


I j gj^j  ^ = Difference  between  the  observed  and  computed  residual 
gravity  anomalies  at  the  i^^  observation  point  for 


iteration. 

g = Observed  residual  gravity  ainomaly  at  the  i''  observation 

i 

point. 

til 

ig  \ = Computed  residual  gravity  anomaly  at  the  i observation 


k i 


point  for  the  k^^  iteration. 


Updated  values  for  the  depths  are  then  computed  using  the  gravity 
differences  as  shown  in  equation  5~28. 


-ki  ^(k-l)i  2 Ti  K p 


(5-28) 


where : 


th  til 

= The  updated  depth  at  the  i observation  point  for  the  k^ 

iteration. 

Tlie  procedure  is  iterated  until  the  gravity  differences  meet  some 
minim\jm  criterion.  Qureshi  and  Mula  (l9Tl)  present  a complete  derivation 
and  discussion  of  the  method  based  on  algorithms  suitable  for  \ise  with 
digital  computers. 
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The  excess  mass,  M,  is  computed  using  the  integral  eqmtion  5-29 


which  is  based  on  Causs'  theorem. 


* 00  • 00 

&y.  (x,y)  dx  dy  = 2 TT  KM 

^ mJX)  j .^oo 


(5-29) 


Equation  5-29  is  entirely  independent  of  shape  and  density. 


In  theory  the  integrations  must  be  carried  over  the  entire  XY  plane. 
However,  the  actual  reinges  of  integrations  are  -X  to  X and  -Y  to  Y 
when  the  origin  is  placed  near  the  center  of  the  residual  gravity  anomaly 
in  the  XY  plane.  The  ranges  in  X and  Y represent  the  limits  of  data 
or  the  distances  at  which  the  residual  gravity  anomalies  become  negligible. 
Equation  5-29  is  rewritten  with  a remainder  term. 


f X f Y 

2 TT  K M = g^  (x,y)  dx  dy  + R {X,Y) 

y_x  j_Y  ^ 


(5-30) 


-X  J-Y 


If  the  center  of  mass  of  M is  at  (x,y,z),  the  remainder  term  is  approx- 


imated by  equation  5-31,  provided  that  |x|  « X and  |y | « Y: 


R(X,Y)  = 27rKM-l+GM  tan-‘  3 


- + y2)l/2 


(5-31) 


Equation  5"32  is  then  derived  by  substituting  equation  5-31  into  equation 


5-30: 


f X r Y 

g„(x,y)  dx  dy 

J _Y  J _Y 


■>-X  •'-Y 
U K tan”^  ~ 


(5-32) 


[X^  + Y^)' 


The  value  for  z is  estimated  by  the  half-width  method  from  equation 
5-2  for  three-dimensional  bodies  (i.e.,  z = 1. 305X1/2)*  Equation  5-5 
is  used  to  compute  z for  two-dimensional  bodies,  z = X1/2  and  equation 


5-32  becomes: 


(x)  dx 


U K tan”  — 


(5-33) 


Simpson’s  Rule  or  some  similar  numerical  quadrature  method  is  used 
to  evaluate  the  integrals  in  equations  5-32  and  5-33  if  the  observed  data 
points  eure  evenly  distributed.  If  the  observed  data  points  are  not  evenly 
distributed,  a template  method  must  be  used  as  described  by  Grant  and 

West  (1965). 

The  inversion  methods  discussed  in  this  section  are  used  to  compute 
first  approximations  to  the  xmknown  parameter  controlling  depth,  size 
and  shape. 

More  complex  techniques  are  beyond  the  scope  and  purpose  of  this 
report.  The  reader  interested  in  such  techniques  is  referred  to  the 
bibliography. 
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6.  CONCLUDING  REMARKS 

The  formulas  and  techniques  presented  in  this  publication  form  the  > 

analytical  basis  of  gravitational  modeling  and  interpretation.  In  appli-  .i 

f 

cation,  several  theoretical  and  practical  aspects  must  be  considered.  < 

i 

The  theoretical  aspects  refer  to  the  inherent  ambiguities  due  to 
the  assumptions  regarding  density  or  mass  distributions  as  well  as  ; 

size  and  shape.  Practical  aspects  to  be  considered  are:  the  amount  i 

j 

and  accuracy  of  the  gravity  data;  the  resolving  power  of  regional-  | 

residual  separation  techniques;  the  amount  of  time  and  mampower  available; 
the  kinds  of  available  digitizing  and  computer  equipment;  the  availa- 
bility and  amoiint  of  the  other  geophysical  data  which  can  be  used  as 
control;  and  most  important,  the  purpose  of  the  investigation. 

Most  systems  of  gravitational  analysis  involve  iterative  combinations 
of  gravitational  modeling  sind  interpretation.  First  approximations 
of  the  unknown  parameters  are  made  using  gravity  and  other  available 
geophysical  data.  The  first  approximations  are  then  used  as  input  in 
appropriate  two-  or  three-dimensional  attraction  equations  and  the  gravi- 
tational effects  are  computed  over  a profile  or  area. 

The  computed  gravitational  effects  are  compared  with  the  observed 
residual  gravity  anomalies  and  adjusted  at  each  computation  point  by 
some  appropriate  mathematical  technique  (i.e.  , analysis  of  Fourier 
components,  Fourier  convolution,  method  of  least  squares  minimization 
or  non-linear  optimization).  The  updated  model,  based  on  the  adjusted 
unknown  parameters , is  then  analyzed  in  terms  of  geological  and  geophysical 
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possibilities.  It  is  often  necessary  and  useful  to  compute  and  analyze 
several  possible  models.  Thus,  the  analytical  techni(iues  of  gravitational 


k 


analysis  provide  generalized  schematic  models.  Detailed  structure 
must  then  be  hypothesized  on  the  basis  of  experience  and  intuitive 
Judgment. 
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